Site Selection for Starbucks¶
Giulia Carella
SDSC20 Workshop - Integrating CARTOframes into Spatial Data Science workflows. October 21 2020
Context¶
Site selection refers to the process of deciding where to open a new or relocate an exiting store/facility by comparing the merits of potential locations. If some attribute of interest is known (e.g. the revenues of existing stores) statistical and machine learning methods can be used to help with the decision process.
What we will do¶
In this demo, we will look at where Starbucks should open new coffee shops in Long Island, NY. Specifically, we will build a model for the revenues of the existing stores as a funcion of sociodemographic data and use this model to predict the expected revenue in each block group.
We will:
- Read and upload to CARTO the revenue data
- Geocode the store adresses (i.e. convert addresses to coordinates)
- Download for Long Island the US Census data by block group
- Enrich each store location with the corresponding block group level Census data
- Build a predictive model for the store revenues
- Publish a visualization with the results
- Discuss how to improve the model results using spatial modelling
Conclusions¶
- Using CARTO we were able to explore and visualize the Starbucks data and enrich them using demographics and socioeconomic data from CARTO Data Observatory
- We built a Generalized Linear model to predict the stores annual revenue from demographic and socioeconomic covariates, which were first transformed to impute missing data and reduce the model dimensionality.
- Predictive checks show that the model is accurate (pseudo R2 ~0.96) and we found that stores with larger annual revenues are located in richer areas and stores with a lower annual revenues are located in areas with a higher workforce density and with a higher concentration of mostly female seniors.
- By using a spatial model we were then able to show how to correct for spatial dependence in the residuals and improve the results.
0. Set up¶
Set up the virtual environment¶
First of all, install virtualenv inside your project folder to secure your project directory to avoid conflict with your other packages.
pip install virtualenv
After installing this run this command one by one inside your root project directory:
virtualenv venv
source venv/bin/activate
Now Your directory is secure and you can install your required packages inside.
pip install -r requirements.txt
To exit the virtual environment type:
deactivate
Import the libraries¶
from requirements import *
from ppca import PPCA
from utils import *
Set CARTO credentials¶
Go to your CARTO dashboard anc click on API KEY to get your key.
carto_username = 'XXXXXX'
carto_API = 'XXXXXX'
set_default_credentials(username=carto_username, api_key=carto_API)
1. Upload Starbucks store data to your CARTO account¶
First we upload the Starbucks data that are stored in your local in the starbucks_long_island.csv file. The data contains the addresses of Starbucks stores, their annual revenue ($) and some store characteristics (e.g. the number of open hours per week).
stores = pd.read_csv('./data/starbucks_long_island.csv')
stores.head(3)
Next, we will plot the distribution of the annual revenues by store.
pandas_hist(stores, 'revenue',bins = 20, textsize = 15, title = 'Annual Revenue ($) by Starbucks store')
2. Geocode locations¶
Since we only know the addresses, we first need to geocode them to extract the corresponding geographical locations, which will be used to assign to each store the Census data of the corresponding block group.
geocoder = Geocoding()
stores = geocoder.geocode(stores, street='addresslines').data
stores.crs = {'init' :'epsg:4326'}
Upload data to CARTO¶
to_carto(stores,'starbucks_long_island_geocoded', if_exists = 'replace')
Visualize the store locations¶
bk = [i*10**5 for i in [4, 7, 10, 13, 16]]
Map(
Layer('starbucks_long_island_geocoded',
style = size_bins_style('revenue',
breaks = bk,
size_range = [1, 30],
color = '#00eaff',
stroke_color = '#006bff'),
legends = size_bins_legend(title='Annual Revenue ($)',
description='STARBUCKS',
footer =''),
),show_info = True,
viewport={'zoom': 10, 'lat': 40.646891, 'lng': -73.869378},
size=(920,300)
)
3. Download demographic and socioeconomic data by block group¶
Next, we will download from CARTO DATA Observatory (www.carto.com/data) the demographic and socioeconomic variables that we will use to build a model for the stores revenues. We will use data from the The American Community Survey (ACS), which is an ongoing survey by the U.S. Census Bureau. The ACS is available at the most granular resolution at the census block group level, with the most recent geography boundaries available for 2015. The data that we will use are from the survey run from 2013 to 2017.
Upload Long Island geometry¶
As we are interested only in the data for the Long Island area, we first upload the geographical boundaries of this area from the manhattan_long_island.geojson file, which is stored locally.
aoi = gpd.read_file('./data/manhattan_long_island.geojson')
Map(
Layer(aoi,
style = basic_style(color = '#eaff00',
stroke_color = '#eaff00',
stroke_width = 2,
opacity = 0.2
),
),show_info = True,
viewport={'zoom': 7.5, 'lat': 40.845419, 'lng': -72.841376},
size=(920,300)
)
Look for the datasets of interest using the Data Discovery methods¶
More documentation here: https://carto.com/developers/cartoframes/guides/Data-discovery/
Let's start by exploring the Data Catalog to see what data categories CARTO can offer for data in the US.
Catalog().country('usa').categories
Then check the available providers and geometries and select the most recent ACS data at the block group level
Catalog().country('usa').category('demographics').datasets.to_dataframe().provider_id.unique()
catalog_usa_demo = Catalog().country('usa').category('demographics').datasets.to_dataframe()
catalog_usa_demo_acs = catalog_usa_demo[catalog_usa_demo.provider_id == 'usa_acs']
catalog_usa_demo_acs.geography_name.unique()
catalog_usa_demo_acs_bk = catalog_usa_demo_acs[catalog_usa_demo_acs.geography_name == 'Census Block Group - United States of America']
catalog_usa_demo_acs_bk.id.to_list()
And explore in greater detail the metadata
Dataset.get('carto-do-public-data.usa_acs.demographics_sociodemographics_usa_blockgroup_2015_5yrs_20132017').to_dict()
Geography.get('carto-do-public-data.carto.geography_usa_blockgroup_2015').to_dict()
Finally get the data and geometries
df_X = Dataset.get('acs_sociodemogr_b758e778')
q = "select * from $dataset$ where ST_Intersects(geom,st_geogfromtext('{}'))".format(aoi.geometry.astype(str)[0])
df_X = df_X.to_dataframe(sql_query=q)
df_X['geometry'] = df_X['geom'].apply(lambda x: str_to_geom(x))
df_X = GeoDataFrame(df_X,geometry = df_X.geometry)
df_X.crs = {'init' :'epsg:4326'}
df_X.head()
Map(
Layer(df_X,
geom_col = 'geometry',
encode_data=False,
style = basic_style(color = '#eaff00',
stroke_width = 1,
opacity = 0.2
),
),show_info = True,
viewport={'zoom': 10, 'lat': 40.646891, 'lng': -73.869378},
size=(920,300)
)
4. Model preparation: normalization and dimensionality reduction¶
First normalize extensive variables by the total population¶
When comparing data for irregular units like census block groups, extra care is needed for extensive variables, i.e. one whose value for a block can be viewed as a sum of sub-block values, as in the case of population. For extensive variables in fact we need first to normalize them, e.g. by dividing by the total area or the total population, depending on the application. Using the metadata available in CARTO Data Observatory we will check which of all the variables can be considered extensive by looking at their default aggregation method: if the aggregation method is classified as SUM, then the variable is normalized by dividing by the total population.
## Get the default aggregation methods
agg_methods_table = read_carto("""
SELECT column_name,
agg_method
FROM variables_public
WHERE dataset_id = 'carto-do-public-data.usa_acs.demographics_sociodemographics_usa_blockgroup_2015_5yrs_20132017'""",
credentials = Credentials(username='do-metadata', api_key='default_public'))
agg_methods_table.dropna(inplace = True)
agg_methods_table.head()
## Compute densities wrt to total population
for i in agg_methods_table.column_name:
if((agg_methods_table[agg_methods_table.column_name==i].agg_method.iloc[0] == 'SUM') & (i!='total_pop')):
df_X[i + '_dens'] = df_X[i].div(df_X['total_pop'])
df_X.drop(i, axis = 1, inplace = True)
## Drop unused columns
drop_cols = df_X.loc[:, df_X.dtypes != np.float64].columns
drop_cols = [i for i in df_X.columns if 'do_' in i or i in drop_cols and i not in ['geoid','geometry']] + ['total_pop']
df_X.drop(drop_cols, axis = 1, inplace = True)
df_X_cols = list(df_X.drop(['geoid','geometry'], axis = 1).columns)
Any missing data?¶
Missing data are common in surveys, as the ACS. Before using the ACS data as covariates we therefore need to check if there are any missing data and to decide how to impute the missing values.
For each variable, the figure below shows the data points that are missing (corresponding to the rows in white). We notice that the amount of missing observations is typically small but depends on the variable.
fig, ax = plt.subplots(figsize=(60,10), nrows =1, ncols = 1)
msno.matrix(df_X[df_X_cols], labels = True, ax = ax)
Are the variables correlated?¶
Next, we will check the pair-wise correlations between the ACS variables. Correlated covariates are problematic for inference in linear (and generalized linear) models as they lead to inflation of the variance of the regression coefficient (https://en.wikipedia.org/wiki/Variance_inflation_factor) but do not necessarly represent a problem for prediction. However, if correlated variables are present, we can take advantage of these correlations to reduce the dimensionality of the model and save computation time.
Looking at the pair-wise correlations of a random sample of 50 variables, we notice large correlations.
draw_corrm(df_X[random.sample(df_X_cols, 50)])
To account for missing values and reduce the model dimensionality, we will transform the data using Priciple Component Analysis¶
To impute the missing data and reduce the model dimensionality we will use a probabilistic formulation of Principal Component Analysis (PCA). PCA is a technique to transform data sets of high dimensional vectors into lower dimensional ones that finds a smaller-dimensional linear representation of data vectors such that the original data could be reconstructed from the compressed representation with the minimum square error.
In its probabilistic formulation, probabilistic PCA (PPCA) the complete data are modelled by a generative latent variable model which iteratively updates the expected complete data log-likelihood and the maximum likelihood estimates of the parameters. In PPCA, the data are incomplete because the principal components are not observed and are treated as latent variables, which also allows to extend the method to when missing data are present.
More information can be found at http://www.jmlr.org/papers/volume11/ilin10a/ilin10a.pdf.
PPCA in a nutshell¶
PCA can also be described as the maximum likelihood solution of a probabilistic latent variable model, which is known as PPCA:
\begin{equation*} Y_{ij} = \mathbf{P}_i \, \mathbf{Z}_j + m_i + \varepsilon_{ij} \quad i = 1, .., d; \, j = 1, .., n \end{equation*}with
\begin{align} p(\mathbf{Z}_j) \sim N(0, \mathbb{1}) \\ p(\varepsilon_{ij}) \sim N(0, \nu) \end{align}Both the principal components $Z$ and the noise $\varepsilon$ are assumed normally distributed. The model can be identified by finding the maximum likelihood (ML) estimate for the model parameters using the Expectation-Maximization (EM) algorithm by minimizing the mean-square error of the observed part of the data. EM is a general framework for learning parameters with incomplete data which iteratively updates the expected complete data log-likelihood and the maximum likelihood estimates of the parameters. In PPCA, the data are incomplete because the principal components, $Z_i$, are not observed and are treated as latent variables. When missing data are present, in the E-step, the expectation of the complete-data log-likelihood is taken with respect to the conditional distribution of the latent variables given the observed variables. In this case, the update EM rules are the following.
E-step: \begin{align} \mathbf{\Sigma}_{\mathbf{Z}_j} = \nu \left(\nu \, \mathbb{1} + \sum_i \mathbf{P}_i \mathbf{P}_i^T \right)^{-1} \\ \overline{\mathbf{Z}}_j = \dfrac{1}{\nu}\mathbf{\Sigma}_{\mathbf{Z}_j} \sum_i \mathbf{P}_i \left(Y_{ij}- m_i \right) \\ m_{i} = \dfrac{1}{n} \sum_j \left(Y_{ij} - \mathbf{P}_i^T \, \overline{\mathbf{Z}}_j \right) \\ \end{align}
M-step: \begin{align} \mathbf{P}_{i} = \left( \sum_j \overline{\mathbf{Z}}_j \overline{\mathbf{Z}}_j ^T + \mathbf{\Sigma}_{\mathbf{Z}_j} \right)^{-1} \sum_j \overline{\mathbf{Z}}_j \, \left(Y_{ij}- m_{ij} \right)\\ \nu = \dfrac{1}{n} \sum_{ij} \left[ \left(Y_{ij} - \mathbf{P}_i^T \, \overline{\mathbf{Z}}_j - m_i \right)^2 + \mathbf{P}_i^T \, \overline{\mathbf{Z}}_j \mathbf{P}_i \, \right] \end{align}
where each row of $\mathbf{P}$ and $\overline{\mathbf{Z}}$ is recomputed based only on those columns of $\overline{\mathbf{Z}}$ which contribute to the reconstruction of the observed values in the corresponding row of the data matrix.
Impute missing values and extract PC scores (linearly decorrelated)¶
X, var_exp = run_ppca(df_X, df_X_cols, min_obs = 0.8)
X.head()
Plot explained variance¶
First, we compute the explained variance as a function of the number of PC scores retained to decide how many PCs we should keep. We decide to keep the PCs up that explain up to 80% of the variance (the black dashed line in the plot below), but this choice is somewhat arbitrary and might vary depending on the application.
plot_pc_var(var_exp, textsize = 15, title = 'PPCA explained variance ratio')
Plot PC correlations¶
To understand the relashionship between the transformed variables (the PC scores) and the original variable, we can plot the 10 variables most highly correlated with each PC, as shown in the bar plot plot which shows the correlation for each of these variables. For example, we see that the first PC, which is that that explains the most of the variance, is positively correlated with the density of ownner occupied housing units but negatively correlated with the density of renter occupied housing units, as also the map below shows.
Note: depending on the random seed on your computer, the signs of the correlations might change.
plot_pc_corr(X, df_X_cols)
Map the first Principal Component¶
map0 = Map(
Layer(X,
encode_data=False,
style = color_bins_style('pc_0',
palette = 'SunsetDark',
stroke_width = 1,
opacity = 0.5),
legends = color_bins_legend(title='First Principal Component Score (PC0)',
description='',
footer ='ACS'),
popup_hover=[popup_element('pc_0', title='First Principal Component Score (PC0)')]
))
map1 = Map(
Layer(X,
encode_data=False,
style = color_bins_style('owner_occupied_housing_units_dens',
palette = 'SunsetDark',
stroke_width = 1,
opacity = 0.5),
legends = color_bins_legend(title='Density of Owner Occupied Housing Units',
description='',
footer ='ACS'),
popup_hover=[popup_element('owner_occupied_housing_units_dens', title='Density of Owner Occupied Housing Units')]
))
map2 = Map(
Layer(X,
encode_data=False,
style = color_bins_style('housing_units_renter_occupied_dens',
palette = 'SunsetDark',
stroke_width = 1,
opacity = 0.5),
legends = color_bins_legend(title='Density of Renter Occupied Housing Units',
description='',
footer ='ACS'),
popup_hover=[popup_element('housing_units_renter_occupied_dens', title='Density of Renter Occupied Housing Units')]
))
Layout([map0,map1,map2], 1,3,viewport={'zoom': 10, 'lat': 40.646891, 'lng': -73.869378})
Select only up the PC up to expalain 80 % of the variance (dimensionality reduction)¶
cum_var_exp = np.cumsum(var_exp)
ncomponents = np.where(cum_var_exp > 0.8)[0][0]
X = X[['geoid','geometry'] + ['pc_' + str(j) for j in range(ncomponents)]]
Finally, Merge the ACS the store data and group the store revenue data by block group¶
stores_enriched = gpd.sjoin(stores[['revenue','the_geom']], X, how='right', op='intersects').drop('index_left', axis = 1)
stores_enriched.head()
5. Train the model using the data for the exisiting stores¶
Having prepared the model covariates, then we will model the annual revenues by store ($\mathbf{y}$) using a Generalized Linear Model (GLM) with the selected PC scores as covariates ($\mathbf{X}$).
\begin{equation*} E\big[\mathbf{y}\vert\mathbf{X}\big] = g^{-1}\big(\sum_j \; X_j * beta_j \big) \end{equation*}where $E$ is the expectation value of the observations $\mathbf{y}$ conditional on the covariates $\mathbf{X}$ and $g$ is a link function. Given the distribution of the annual revenues, we will use a Gamma likelihood with a logarithmic link function.
df1 = stores_enriched[~np.isnan(stores_enriched['revenue'])]
lin_scale = 10**5
pc_cols = [i for i in df1.columns if 'pc_' in i]
y = df1['revenue'].values/lin_scale
X = df1[pc_cols]
X0 = np.ones((X.shape[0],1))
X = np.hstack((X,X0))
sm_fit = sm.GLM(y,X, family=sm.families.Gamma(link=sm.families.links.log)).fit()
df1['y_hat'] = sm_fit.predict(X)*lin_scale
Looking at the model results we notice a significant negative correlation with the first PC (which is lower in richer areas) and a significant positive correlation with the second and third PCs which is negatively correlated with a higher workforce density and a higher density of female seniors.
Note: the individual correlations with the PC and therefore with the original variables, might depend on the random seed on your computer but will be consistent overall.
print(sm_fit.summary())
To asses the model accuracy, we then plot the observed vs. the predicted values and compute the pseudo-$R^2$ score which is defined as
\begin{equation*} 1 - \dfrac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \bar{y})^2} \end{equation*}As shown by the plot, below which indicates that most of the data points lie on the 1-1 line and a high value of the pseudo-$R^2$ score, we can conclude that our model is accurate in the in-sample predictions.
pandas_scatter(df1,'revenue','y_hat', xlab = 'Observed', ylab = 'Predicted', pseudo_R2 = True)
While GLM model typically do not overfit, we can also check the pseudo-$R^2$ score of the training and testing errors.
X_train, X_test, y_train, y_test, idx_train, idx_test \
= train_test_split(X, y, np.arange(X.shape[0]), test_size=0.2, random_state=42)
fitglm = sm.GLM(y_train,X_train, family=sm.families.Gamma(link=sm.families.links.log)).fit()
df1_train = df1.iloc[idx_train]
df1_test = df1.iloc[idx_test]
print('Training pseudo-r2 score: ' + str(pseudo_r2score(y_train,fitglm.predict(X_train))))
print('Training pseudo-r2 score: ' + str(pseudo_r2score(y_test,fitglm.predict(X_test))))
6. Get the predictions by Census block group¶
Finally, we can use this model to predict the annual revenue ($) for each block group in Long Island.
df2 = stores_enriched.copy()
y = df2['revenue'].values/lin_scale
X = df2[pc_cols]
X0 = np.ones((X.shape[0],1))
X = np.hstack((X,X0))
df2['revenue_pred'] = sm_fit.predict(X)*lin_scale
7. Visualize and share¶
With CARTOframes you can then share your visualizations by publishing your map using a url.
Upload data to CARTO¶
to_carto(df2,'starbucks_long_island_predictions', if_exists = 'replace')
Visualize the predictions¶
bk = [i*10**5 for i in [4, 7, 10, 13, 16, 20, 30, 40]]
pubmap = Map(
[Layer('starbucks_long_island_predictions',
style = color_bins_style('revenue_pred',
breaks = bk,
palette = 'SunsetDark',
stroke_width = 1,
opacity = 0.5),
legends = color_bins_legend(title='Annual Predicted Revenue ($)',
description='',
footer =''),
popup_hover=[popup_element('revenue_pred', title='Annual Predicted Revenue ($)')],
widgets=[histogram_widget('revenue_pred',
title='Annual Predicted Revenue ($)',
description='Select a range of values to filter',
buckets=15)]),
Layer('starbucks_long_island_geocoded',
style = size_bins_style('revenue',
breaks = bk,
size_range = [1, 30],
color = '#00eaff',
stroke_color = '#006bff'),
legends = size_bins_legend(title='Annual Revenue ($)',
description='STARBUCKS',
footer ='2019'))
],show_info = True,
viewport={'zoom': 9.4, 'lat': 40.673790, 'lng': -73.992685},
size=(920,300)
)
pubmap
Share your results¶
pubmap.publish("starbucks_long_island",
password=None,
if_exists='replace')
5. Can we improve our model by accounting for any spatial dependence in the residuals?¶
Having assessed the performance of the non-spatial GLM model, we can check if the model residuals ${y_i-\hat{y_i}}$ show any residual spatial dependence, indicating that better results might be obtained using a model that explicitly takes the property that an observation is more correlated with an observation collected at a neighboring location than with another observation that is collected from farther away.
Including a spatially-structured random effects model that incorporates such spatial dependency is beneficial when the variables used as covariates in the model may not be sufficient to explain the variability of the observations and the measurements, given the predictors, are not independent. Including a spatially-structured random effect in the model helps mitigating the situation.
How to go spatial for spatially continuous processes?¶
A spatial stochastic process can be decomposed as:
\begin{equation*} y(\mathbf{s}) = \mu(\mathbf{s}) + w(\mathbf{s}) + \varepsilon \end{equation*}where $\mu(\mathbf{s})$ is the mean function (deterministic and smooth), $\varepsilon$ is an independent process and
\begin{equation*} w(\mathbf{s}) \sim GP\left(0, C(\mathbf{s}, \mathbf{s'}, \mathbf{\theta})\right) \end{equation*}where the covariance $C(\mathbf{s}, \mathbf{s'}, \mathbf{\theta})$ typically depends on the Euclidean length of the spatial separation vector only $|\mathbf{h}| = |\mathbf{s} - \mathbf{s'}|$, and potentially on some hyperparameters $\mathbf{\theta}$. For example the exponential covariance function is defined as:
- Exponential covariance:
where the empirically derived definition for the range is $\rho = \dfrac{\sqrt{8 \nu}}{\kappa}$.
with the covariance between two points tending to zero as the points become further separated in space.
Fitting spatial models, comes with some caveats:
- Fitting the above model requires calculating determinants and inverses of large matrices: dense matrices scales cubically with matrix size, given by the number of locations where the process is observed
- Computational complexity grows with matrix size matrix complexity
- Memory limitations also create problems
Variograms¶
Variograms represent a useful tool to check for spatial dependence in the residuals. The variogram is defined as:
\begin{equation*} Var\left[y(\mathbf{s + h}) - y(\mathbf{s})\right]^2 =: 2\gamma(\mathbf{h}) = 2 \left(C \left(0\right) - C \left(|\mathbf{h}|\right)\right)\end{equation*}First, we will construct an empirical variogram by grouping the pairs of observations in bins based on their distance (if isotropy is assumed) and averaging the squared differences from the values for all pairs. Given a semivariogram model, we can then fit this model to the empirical semivariogram and estimate the model parameters.
First, fit a variogram model on the residuals¶
coords_x = df1['geometry'].apply(lambda x: x.centroid.x)
coords_y = df1['geometry'].apply(lambda x: x.centroid.y)
coordinates = pd.concat([coords_x, coords_y], axis=1).values
values = (df1['y_hat']-df1['revenue'])/lin_scale
vg = Variogram(coordinates=coordinates,
values=values,
normalize=False,
use_nugget=True,
model='exponential',
maxlag="mean",
n_lags=15)
ax = Variogram_plot(vg, fig_title='Residual variogram')
The variogram models parameters are:
print(vg)
Then, let's fit a Gaussian Process GLM model¶
We will now extend the GLM model to account for the residual spatial dependence, by fitting this model:
\begin{equation*} E\big[\mathbf{y}\vert\mathbf{X}\big] = g^{-1}\,\Big(\sum_j \; X_j * beta_j + u(\mathbf{s}_i) + \varepsilon_i \Big) \end{equation*}where $u(\mathbf{s}_i)$ represents the spatially structured residuals modelled as a Gaussian Process with an exponential covariance function and $\varepsilon_i$ ad IID random term.
We will fit this model in a Bayesian framework using the programming language Stan (https://mc-stan.org). However, to reduce the model complexity we will fix the covariance parameters and use those derived from the variogram analysis (in a fully Bayesian approach these would be estimated).
The stan model is saved locally in _stan/GP_gamma_modelexact.stan.
y = stores_enriched['revenue'].values/lin_scale
X = stores_enriched[pc_cols]
X0 = np.ones((X.shape[0],1))
X = np.hstack((X,X0))
coords_x = stores_enriched['geometry'].apply(lambda x: x.centroid.y)
coords_y = stores_enriched['geometry'].apply(lambda x: x.centroid.y)
coords = np.vstack((coords_x.values,coords_y.values)).T
y1 = y[~np.isnan(y)]
X1 = X[~np.isnan(y)]
coords1_x = df1['geometry'].apply(lambda x: x.centroid.y)
coords1_y = df1['geometry'].apply(lambda x: x.centroid.y)
coords1 = np.vstack((coords1_x.values,coords1_y.values)).T
pho, sigma, tau = vg.parameters
## Define the data dictionary the model
data = {'N': X1.shape[0],'K': X1.shape[1],
'X': X1, 'y': y1, 'coords':coords1,
'pho':pho, 'sigma':sigma, 'tau':tau}
## Compile the model
#stanm = pystan.StanModel(file='./stan/GP_gamma_model_exact.stan')
## Train the model and generate samples
#stan_fit = stanm.sampling(data=data, iter=1000, chains=3,
# warmup=500, thin=1, seed=101)
## Save the model
#with gzip.open("./stan/GP_gamma_model_exact.gz", "wb") as f:
# pickle.dump({'model' : stanm, 'fit' : stan_fit}, f, protocol=-1)
## Load the model
with gzip.open("./stan/GP_gamma_model_exact.gz", "rb") as f:
data_dict = pickle.load(f)
stan_fit = data_dict['fit']
Extract the posterior mean and standard deviations
df1['y_hat_m'] = stan_fit.extract()['y_rep'].mean(axis=0)*lin_scale
df1['y_hat_se'] = stan_fit.extract()['y_rep'].std(axis=0)*lin_scale
pandas_scatter(df1,'revenue','y_hat_m', xlab = 'Observed', ylab = 'Predicted', pseudo_R2 = False)
_ = plt.plot([df1.revenue,df1.revenue],
[df1['y_hat_m']-df1['y_hat_se'], df1['y_hat_m']+df1['y_hat_se']], 'b-', linewidth = 0.3,linestyle = '--')
plt.xlim((df1['y_hat_m']-df1['y_hat_se']).min(),(df1['y_hat_m']+df1['y_hat_se']).max())
Finally check that we have removed any residual spatial autocorrelation¶
We fit now the variogram mode to the model residuals, using the posterior mean.
Compared to the previous variogram plot, we can see that by adding a spatially structured random effect we have consistently reduced the spatial dependence in the residuals (check the scale of the y-axis).
## Fit variogram model
values = (df1['y_hat_m']-df1['revenue'])/lin_scale
coordinates = coords1
vg = Variogram(coordinates=coordinates,
values=values,
normalize=False,
use_nugget=True,
model='exponential',
maxlag="mean",
n_lags=15)
ax = Variogram_plot(vg, fig_title='Residual variogram')
Next steps:¶
Clearly using the above model we have showed that we have removed the spatial dependence in the residuals. The next step would then be to use this model to make the predictions for the annual revenues for each census block group. However, because the computation time for GPs scales cubically with the data locations, we need some approximations which goes beyond the scope of this workshop. More information can be found here:
- Using Stan: https://arxiv.org/abs/2004.11408
- Using the Integrated Nested Laplace Approximation (INLA) with the stochastic partial differential equation approach (SPDE): http://www.r-inla.org, https://becarioprecario.bitbucket.io/spde-gitbook/, https://rss.onlinelibrary.wiley.com/doi/full/10.1111/j.1467-9868.2011.00777.x,